In [1]:
%load_ext load_style
%load_style talk.css
This notebook is inspired by the NCL SPI example.
The Standardized Precipitation Index (SPI) is a probability index that gives a better representation of abnormal wetness and dryness than the Palmer Severe Drought Index (PSDI). The World Meteorological Organization (WMO) recommends, that all national meteorological and hydrological services should use the SPI for monitoring of dry spells. Some advantages of the SPI:
A shortcoming of the SPI, as noted by Trenbert et al (2014):
In this notebook, SPI is obtained by fitting a gamma distribution to monthly GPCP precipitation data from 1979 to 2010. The data can be downloaded from https://www.ncl.ucar.edu/Applications/Data/.
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt # to generate plots
from mpl_toolkits.basemap import Basemap # plot on map projections
import datetime
from netCDF4 import Dataset # http://unidata.github.io/netcdf4-python/
from netCDF4 import netcdftime
from netcdftime import utime
from dim_spi_n import * # private lib
from timeit import default_timer as timer
import warnings
warnings.simplefilter('ignore')
In [3]:
infile = r'data/V22_GPCP.1979-2010.nc'
fh = Dataset(infile, mode='r') # file handle, open in read only mode
fh.set_auto_mask(False)
lons = fh.variables['lon'][:]
lats = fh.variables['lat'][:]
nctime = fh.variables['time'][:]
t_unit = fh.variables['time'].units
pr = fh.variables['PREC'][:]
try :
t_cal = fh.variables['time'].calendar
except AttributeError : # Attribute doesn't exist
t_cal = u"gregorian" # or standard
fh.close() # close the file
undef = -99999.0
pr[pr==undef] = np.nan
pr = pr.astype(np.float64)
nt,nlat,nlon = pr.shape
ngrd = nlat*nlon
In [4]:
utime = netcdftime.utime(t_unit, calendar = t_cal)
datevar = utime.num2date(nctime)
datevar[0:5]
Out[4]:
In [5]:
pr_grd = pr.reshape((nt,ngrd), order='F')
spi_grd = np.zeros(pr_grd.shape)
spi_grd[:,:] = np.nan
nrun = 24
s = timer()
for igrd in np.arange(ngrd):
one_pr = pr_grd[:,igrd]
if (isinstance(one_pr, np.ma.MaskedArray)) and one_pr.mask.all():
print(igrd)
continue
else:
spi_grd[:,igrd] = dim_spi_n(one_pr, nrun)
spi = spi_grd.reshape((nt,nlat,nlon), order='F')
e = timer()
print(e - s)
In [6]:
m = Basemap(projection='cyl', llcrnrlon=min(lons), llcrnrlat=min(lats),
urcrnrlon=max(lons), urcrnrlat=max(lats))
x, y = m(*np.meshgrid(lons, lats))
clevs = np.linspace(-3.0, 3.0, 21)
fig = plt.figure(figsize=(15,12))
# Plot the first one
ax = fig.add_subplot(211)
idx = 258
cs = m.contourf(x, y,spi[idx,:,:], clevs, cmap=plt.cm.RdBu)
m.drawcoastlines()
cb = m.colorbar(cs)
plt.title('SPI-' + str(nrun) + ' In '+ datetime.date.strftime(datevar[idx], "%m/%Y"), fontsize=16)
# plot the second one
ax = fig.add_subplot(212)
idx = -1
cs = m.contourf(x, y,spi[idx,:,:], clevs, cmap=plt.cm.RdBu)
m.drawcoastlines()
cb = m.colorbar(cs)
plt.title('SPI-' + str(nrun) + ' In '+ datetime.date.strftime(datevar[idx], "%m/%Y"), fontsize=16)
Out[6]:
McKee, T.B., N.J. Doesken, and J. Kleist, 1993. The relationship of drought frequency and duration ot time scales. Eighth Conference on Applied Climatology, American Meteorological Society Jan 17-23, 1993, Anaheim CA, pp. 179-186.
McKee, T.B., N.J. Doesken, and J. Kleist, 1995. Drought monitoring with multiple time scales. Ninth Conference on Applied Climatology, American Meteorological Society Jan 15-20, 1995, Dallas TX, pp. 233-236.
Guttman, N.B., 1998. Comparing the Palmer Drought Index and the Standardized Precipitation Index. Journal of the American Water Resources Association, 34(1), 113-121.
Guttman, N.B., 1999. Accepting the Standardized Precipitation Index: A calculation algorithm. Journal of the American Water Resources Association, 35(2), 311-322.
Trenberth et al (2014) Global warming and changes in drought Nature Climate Change 4, 17-22; doi:10.1038/nclimate2067
In [ ]: